function IdSet=OneSide3(PYX_cond, u0grid, u1grid, u2grid,...
                        sgb, ...
                        equalorder, ...
                        Ugridreduced, sgreduced,...
                        id_Ugridreduced,...
                        s_id_gridreduced, sU_id_gridreduced,...
                        n_U_sets, Tcoord, indices_pairs,...
                        numeqconstraints,numzeros,coordrowseq,d3, e5, f5, g3, h5, i5, l3,m5, n5, fillAeq)
%% Complete equality constraints for parts invariant across parameter values
beq=[PYX_cond(1)-1;... %observational equivalence
      PYX_cond(2);...
      PYX_cond(3);...
      zeros(numzeros,1);... %zeros
      1;... %ones       
      ones((3+sgb*3)*3,1); ... %symmetry
      1/2*ones(3,1);...
      zeros((2*7)+(2*2*sgb)*3,1)]; %identical  

%% Loop over parameter values 
exists=zeros(sgreduced,1);
eps=10^(-6);
for g=1:sgreduced
% Set indices to keep track of equal elements     
idU1=id_Ugridreduced{1}(g,:);
idU2=id_Ugridreduced{2}(g,:);
idU3=id_Ugridreduced{3}(g,:);
idU1_t=id_Ugridreduced{1}(g,:).';
idU2_t=id_Ugridreduced{2}(g,:).';
idU3_t=id_Ugridreduced{3}(g,:).';
% Relevant sizes of the vector of unknowns
s1=s_id_gridreduced(g,1);
s2=s_id_gridreduced(g,2);
%s3=s_id_gridreduced(g,3);
sU=sU_id_gridreduced(g); 

%%%%%%%%%%%%%%%%%%%%%% Equality constraints %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Coordcolumnseq
id_01=[changecoord(s1,s2,idU1((1:9+3*2*sgb)),idU2(9),idU3(1))...  
       changecoord(s1,s2,idU1(9),idU2((1:8)),idU3(1))...
       changecoord(s1,s2,idU1(9),idU2(10:9+3*2*sgb),idU3(1))...
       changecoord(s1,s2,idU1(repmat((1:1:9+3*2*sgb),1,7+2*sgb*3)), idU2(9), idU3(kron((3:1:9+3*2*sgb), ones(1,9+3*2*sgb))))...
       changecoord(s1,s2,idU1(9),idU2(repmat((1:8),1,7+2*sgb*3)),idU3(kron((3:1:9+3*2*sgb), ones(1,8))))...
       changecoord(s1,s2,idU1(9),idU2(repmat((10:9+3*2*sgb),1,7+2*sgb*3)),idU3(kron((3:1:9+3*2*sgb), ones(1,2*sgb*3))))...
       changecoord(s1,s2,idU1(repmat((1:1:9+3*2*sgb),1,9+3*2*sgb)),idU2(kron((1:1:9+3*2*sgb), ones(1,9+3*2*sgb))),idU3(2))...
       changecoord(s1,s2,idU1(8),idU2(8),idU3(1))]; %zeros and ones


coordcolumnseq= [changecoord(s1,s2,idU1(1),idU2(8),idU3(3)) changecoord(s1,s2,idU1(8),idU2(8),idU3(3)) changecoord(s1,s2,idU1(1),idU2(8),idU3(1)) ... %observational equivalence
                 changecoord(s1,s2,idU1(8),idU2(8),idU3(3)) changecoord(s1,s2,idU1(8),idU2(1),idU3(3))...
                 changecoord(s1,s2,idU1(1),idU2(1),idU3(1)) ...
                 id_01 ...
                 changecoord(s1,s2,idU1(1),idU2(8),idU3(1)) changecoord(s1,s2,idU1(2),idU2(8),idU3(1)) ... %symmetry
                 changecoord(s1,s2,idU1(3),idU2(8),idU3(1)) changecoord(s1,s2,idU1(4),idU2(8),idU3(1)) ...
                 changecoord(s1,s2,idU1(5),idU2(8),idU3(1)) changecoord(s1,s2,idU1(6),idU2(8),idU3(1)) ...
                 changecoord(s1,s2,idU1((10:1:9+3*2*sgb)),idU2(8),idU3(1))...
                 changecoord(s1,s2,idU1(8),idU2(1),idU3(1)) changecoord(s1,s2,idU1(8),idU2(2),idU3(1)) ...
                 changecoord(s1,s2,idU1(8),idU2(3),idU3(1)) changecoord(s1,s2,idU1(8),idU2(4),idU3(1)) ...
                 changecoord(s1,s2,idU1(8),idU2(5),idU3(1)) changecoord(s1,s2,idU1(8),idU2(6),idU3(1)) ...
                 changecoord(s1,s2,idU1(8),idU2((10:1:9+3*2*sgb)),idU3(1))...
                 changecoord(s1,s2,idU1(8),idU2(8),idU3(3)) changecoord(s1,s2,idU1(8),idU2(8),idU3(4)) ...
                 changecoord(s1,s2,idU1(8),idU2(8),idU3(5)) changecoord(s1,s2,idU1(8),idU2(8),idU3(6)) ...
                 changecoord(s1,s2,idU1(8),idU2(8),idU3(7)) changecoord(s1,s2,idU1(8),idU2(8),idU3(8)) ...
                 changecoord(s1,s2,idU1(8),idU2(8),idU3((10:1:9+3*2*sgb)))...
                 changecoord(s1,s2,idU1(7),idU2(8),idU3(1)) ... %1/2
                 changecoord(s1,s2,idU1(8),idU2(7),idU3(1)) ...
                 changecoord(s1,s2,idU1(8),idU2(8),idU3(9)) ...
                 changecoord(s1,s2,idU1(1),idU2(8),idU3(1)) changecoord(s1,s2,idU1(8),idU2(3),idU3(1)) ... %identical
                 changecoord(s1,s2,idU1(1),idU2(8),idU3(1)) changecoord(s1,s2,idU1(8),idU2(8),idU3(5)) ...
                 changecoord(s1,s2,idU1(2),idU2(8),idU3(1)) changecoord(s1,s2,idU1(8),idU2(4),idU3(1)) ...
                 changecoord(s1,s2,idU1(2),idU2(8),idU3(1)) changecoord(s1,s2,idU1(8),idU2(8),idU3(6)) ...
                 changecoord(s1,s2,idU1(3),idU2(8),idU3(1)) changecoord(s1,s2,idU1(8),idU2(1),idU3(1)) ...
                 changecoord(s1,s2,idU1(3),idU2(8),idU3(1)) changecoord(s1,s2,idU1(8),idU2(8),idU3(7)) ...
                 changecoord(s1,s2,idU1(4),idU2(8),idU3(1)) changecoord(s1,s2,idU1(8),idU2(2),idU3(1)) ...
                 changecoord(s1,s2,idU1(4),idU2(8),idU3(1)) changecoord(s1,s2,idU1(8),idU2(8),idU3(8)) ...
                 changecoord(s1,s2,idU1(5),idU2(8),idU3(1)) changecoord(s1,s2,idU1(8),idU2(5),idU3(1)) ...
                 changecoord(s1,s2,idU1(5),idU2(8),idU3(1)) changecoord(s1,s2,idU1(8),idU2(8),idU3(3)) ...
                 changecoord(s1,s2,idU1(6),idU2(8),idU3(1)) changecoord(s1,s2,idU1(8),idU2(6),idU3(1)) ...
                 changecoord(s1,s2,idU1(6),idU2(8),idU3(1)) changecoord(s1,s2,idU1(8),idU2(8),idU3(4)) ...
                 changecoord(s1,s2,idU1(7),idU2(8),idU3(1)) changecoord(s1,s2,idU1(8),idU2(7),idU3(1)) ...
                 changecoord(s1,s2,idU1(7),idU2(8),idU3(1)) changecoord(s1,s2,idU1(8),idU2(8),idU3(9)) ...
                 changecoord(s1,s2,idU1(d3(:).'), idU2(e5(:).'), idU3(f5(:).')) ...
                 changecoord(s1,s2,idU1(g3(:).'), idU2(h5(:).'), idU3(i5(:).')) ...
                 changecoord(s1,s2,idU1(l3(:).'), idU2(m5(:).'), idU3(n5(:).')) ];           

Aeq=sparse(coordrowseq, coordcolumnseq, fillAeq, numeqconstraints, sU);   %[#equalityconstraints]x[sU]

%%%%%%%%%%%%%%%%%%%%%%  Inequality constraints %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Repeat the way we constructed Tcoord_specific but now using the parameter values under consideration
vectors = cellfun(@(x) {x(g,:)}, Ugridreduced); %cell 1xn_U_sets
%For j=1,...,n_U_sets, it picks the g-th row of Ugridreduced{j}
T_temp = cell(1,n_U_sets);
[T_temp{:}] = ndgrid(vectors{:}); 
T_temp = cat(n_U_sets+1, T_temp{:});
T = reshape(T_temp,[],n_U_sets); %matrix of dimension as Tcoord reporting the elements of all the possible n_U_sets-tuples from the U-sets
D=[T(indices_pairs(:,1),:) T(indices_pairs(:,2),:) ...
   Tcoord(indices_pairs(:,1),:) Tcoord(indices_pairs(:,2),:)]; %n_U_sets-tuples n_U_sets-tuples coordinates coordinates
    
%Save in D1 the rows of D_temp where the first triplet is <= then the second
D1_all=D(sum(D(:,1:3)<=D(:,4:6),2)==3, :);
D1=D1_all(:, 7:end); %only coordinates
sd1=size(D1,1);
%Remember <= means: (<,<,<) or (<,=,=) or (<,<,=) or (=,=,=) [any order]
%Identify zero boxes
D1_values= D1_all(:, 1:6);
C1=(round(D1_values(:,1),8)>=round(D1_values(:,5)+D1_values(:,6),8) | round(D1_values(:,4),8)<=round(D1_values(:,2)+D1_values(:,3),8));


%Save in D2 the rows of D_temp where the first triplet is >= (excluding = = =) then the second
D2_all=D((sum(D(:,1:3)>D(:,4:6),2)==3) |... %>,>,>
       (sum(D(:,1:3)>D(:,4:6),2)==2 & sum(D(:,1:3)==D(:,4:6),2)==1) | ... %>,>,=
       (sum(D(:,1:3)>D(:,4:6),2)==1 & sum(D(:,1:3)==D(:,4:6),2)==2), :); %>,=,= 
D2=D2_all(:, 7:end); %only coordinates   
sd2=size(D2,1);
%Identify zero boxes
D2_values=[D2_all(:, 4:6) D2_all(:, 1:3)] ; 
C2=(round(D2_values(:,1),8)>=round(D2_values(:,5)+D2_values(:,6),8) | round(D2_values(:,4),8)<=round(D2_values(:,2)+D2_values(:,3),8));


coordrowsineq=[kron((1:1:sd1), ones(1,8))  kron((sd1+1:1:sd1+sd2), ones(1,8))]; %1x[(sd1+sd2)*8]
%8 comes from the number of all possible vertices generated by each pair of triplets
coordcolumns1=[changecoord(s1,s2,idU1_t(D1(:,1)),idU2_t(D1(:,2)),idU3_t(D1(:,3))) ... 
                   changecoord(s1,s2,idU1_t(D1(:,4)),idU2_t(D1(:,2)),idU3_t(D1(:,3))) ...
                   changecoord(s1,s2,idU1_t(D1(:,1)),idU2_t(D1(:,5)),idU3_t(D1(:,3))) ...
                   changecoord(s1,s2,idU1_t(D1(:,4)),idU2_t(D1(:,5)),idU3_t(D1(:,3))) ...
                   changecoord(s1,s2,idU1_t(D1(:,1)),idU2_t(D1(:,2)),idU3_t(D1(:,6))) ...
                   changecoord(s1,s2,idU1_t(D1(:,4)),idU2_t(D1(:,2)),idU3_t(D1(:,6))) ...
                   changecoord(s1,s2,idU1_t(D1(:,1)),idU2_t(D1(:,5)),idU3_t(D1(:,6))) ...
                   changecoord(s1,s2,idU1_t(D1(:,4)),idU2_t(D1(:,5)),idU3_t(D1(:,6)))]; %sd1x8
coordcolumns2=[changecoord(s1,s2,idU1_t(D2(:,1)),idU2_t(D2(:,2)),idU3_t(D2(:,3))) ... 
                   changecoord(s1,s2,idU1_t(D2(:,4)),idU2_t(D2(:,2)),idU3_t(D2(:,3))) ...
                   changecoord(s1,s2,idU1_t(D2(:,1)),idU2_t(D2(:,5)),idU3_t(D2(:,3))) ...
                   changecoord(s1,s2,idU1_t(D2(:,4)),idU2_t(D2(:,5)),idU3_t(D2(:,3))) ...
                   changecoord(s1,s2,idU1_t(D2(:,1)),idU2_t(D2(:,2)),idU3_t(D2(:,6))) ...
                   changecoord(s1,s2,idU1_t(D2(:,4)),idU2_t(D2(:,2)),idU3_t(D2(:,6))) ...
                   changecoord(s1,s2,idU1_t(D2(:,1)),idU2_t(D2(:,5)),idU3_t(D2(:,6))) ...
                   changecoord(s1,s2,idU1_t(D2(:,4)),idU2_t(D2(:,5)),idU3_t(D2(:,6)))]; %sd2x8
          
coordcolumnsineq=[reshape(coordcolumns1.', 1, 8*sd1) reshape(coordcolumns2.', 1, 8*sd2)];   %transform coordcolumns1,2 in a row vector by reshaping row-wise  
%1x[(sd1+sd2)*8]


fillAineq=[repmat([1 -1 -1 1 -1 1 1 -1], 1, sd1) repmat([-1 1 1 -1 1 -1 -1 1], 1, sd2)]; %remember: signs are flipped because inequalities are <= 
%1x[(sd1+sd2)*8]

bineq=zeros(sd1+sd2, 1);


Aineq=sparse(coordrowsineq, coordcolumnsineq, fillAineq, sd1+sd2, sU);   %[#inequalityconstraints]x[sU]

lower_bineq=-Inf(size(Aineq,1),1);
lower_bineq([C1; C2])=0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Run LP %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


lower=zeros(sU,1)+eps;
upper=ones(sU,1)-eps;
lower(id_01)=0;
upper(id_01)=1;
% MOSEK
param_MOSEK.MSK_IPAR_LOG = 0; 
prob.blx=lower; %lower bound unknowns
prob.ulx=upper; %upper bound unknowns
prob.c=zeros(sU,1); %objective function
prob.a=[Aeq; Aineq];
prob.blc=[beq; lower_bineq];
prob.buc=[beq; bineq];
[~,result]     = mosekopt('minimize echo(0)',prob, param_MOSEK);
if strcmp(result.sol.itr.prosta, 'PRIMAL_AND_DUAL_FEASIBLE') %if the primal is feasible
   exists(g)=1; 
elseif strcmp(result.sol.itr.prosta, 'UNKNOWN') && strcmp(result.sol.bas.prosta, 'PRIMAL_AND_DUAL_FEASIBLE')
   exists(g)=1; 
end



end



%Enlarge results to the original grid of parameter vectors
exists=exists(equalorder);
IdSet = [u0grid(logical(exists),:)  u1grid(logical(exists),:)  u2grid(logical(exists),:)];

